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We discuss the use of Langevin molecular dynamics in the investigation of the non-equilibrium properties 
of disordered vortex matter. Our special focus is set on values of system parameters that are realistic 
for disordered high-Tc superconductors such as YBCO. Using a discretized elastic line model, we study 
different aspects of vortices far from thermal equilibrium. On the one hand we investigate steady-state 
properties of driven magnetic flux lines in a disordered environment, namely the current-voltage charac¬ 
teristics, the gyration radius, and the pinning time statistics. On the other hand we study the complex 
relaxation processes and glassy-like dynamics that emerge in type-11 superconductors due to the intricate 
competition between the long-range vortex-vortex repulsion and flux pinning due to randomly placed 
point defects. To this end we consider different types of sudden perturbations: temperature, magnetic 
field, and external current quenches. 

Keywords: Langevin molecular dynamics; disordered vortex matter; type-ll superconductors; pinning 
statistics; relaxation processes 


1. Introduction 

Langevin equations are stochastic differential equations that have seen widespread applications in 
many areas of physics [1]. They form a sound effective description in situations where a strict sep¬ 
aration of time scales prevails, i.e., where one has both fast and slow degrees of freedom. Langevin 
equations then provide equations of motion for the relevant slow modes, whereas the effects of the 
fast degrees (e.g., the environment) are incorporated via stochastic forcing. 

Since the initial use of a Langevin equation for the description of a Brownian particle [2], many 
other problems have been tackled using stochastic differential equations. Obvious problems can 
be found in biological settings where large ’particles’ are moving in an aqueous environment. 
But there are also important problems in condensed matter and material physics that can be 
successfully addressed using Langevin equations. Well-known examples include domain walls in 
magnetic systems, fluctuating interfaces, or flux lines in type-II superconductors. 

In this paper we discuss the use of Langevin equations for the study of interacting vortices in a 
disordered environment. Our main focus is on high-Tc superconductors [3] where glassy states are 
realized at low temperatures due to the competition of long-range vortex-vortex interaction and 
short-range defect pinning. More specifically, we will discuss in some detail the properties of these 
systems when they are forced far away from thermal equilibrium. By employing the methods of 
Langevin molecular dynamics, we numerically solve the set of coupled stochastic differential equa- 
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tions describing our systems. In particular, our model systems are well suited for the application of 
this method due to the separation into many slow degrees of freedom embedded in an environment 
dominated by fast flluctuating thermal degrees of freedom. 

In principle, a system can be prepared ’far from equilibrium’ in two different ways. In the first 
case the system is subjected to an external drive that yields a non-equilibrium steady state. Non¬ 
equilibrium steady states have been a focus of many studies in statistical physics that aim at finding 
a comprehensive theoretical framework allowing to determine the stationary probability density 
and probability current distribution. Whereas some progress has been made in certain instances, 
a comprehensive theory remains elusive. Non-equilibrium steady states are also of interest in the 
context of vortex matter as they emerge naturally when applying a driving force to the system. 

Another way to induce a system out of equilibrium is to force it away from a stationary state 
through a sudden change of at least one of the system parameters [4]. This sudden change, often 
called quench, can be realized in different ways. Very common is a temperature quench, but other 
types of quenches (magnetic field or vortex density as well as external current or drive quenches) 
are also possible and will be discussed in this paper. Once the system has been brought out of the 
initial steady state, it will relax towards a new stationary state. This relaxation process can be very 
long-lasting in systems with slow or glassy-like dynamics. In some instances this relaxation process 
manifests itself through characteristic dynamical scaling behavior [5]. In other cases relaxation is 
a very involved process, due to the competition between different physical interactions or other 
effects, and time-dependent quantities may exhibit a very rich and complex behavior. The intricate 
slow relaxation processes of interacting magnetic vortices in disordered type-II superconductors 
deep inside the low-temperature Bragg glass phase are of this second type. 

This paper is meant to provide a didactic introduction to the investigation of non-equilibrium 
properties of flux lines in a disordered environment using Langevin molecular dynamics. After in¬ 
troducing the basics of Langevin dynamics in the next section, we focus in section 3 on the study of 
interacting vortex matter using Langevin molecular dynamics. We first recall the effective descrip¬ 
tion through an elastic line model before describing the implementation of Langevin molecular 
dynamics for this system. Section 4 presents a range of results, both for driven vortices in the 
steady state as well as for non-equilibrium relaxation processes that result from a quench of one of 
the external system control parameters. We conclude in section 5. 


2. Langevin dynamics 

Before we employ Langevin molecular dynamics for the investigation of the non-equilibrium prop¬ 
erties of vortices in a disordered environment in the next section, it is useful to first illustrate the 
concept of a Langevin equation through simpler examples. For that reason we will briefly recall 
the case of a Brownian particle, followed be a general discussion of the Langevin equation of a sys¬ 
tem characterized by a Hamiltonian and the associated Boltzmann-Gibbs equilibrium probability 
distribution. 


2.1. Brownian particle 

The paradigmatic example of a Langevin equation is of course the stochastic equation that describes 
a particle immersed in a fluid. This particle is continuously exposed to collisions with the molecules 
of the fluid. Assuming that there is no net flow, the molecules hit the particle with the same 
rate from all sides, i.e., on average the force on the particle due to the fluid molecules is zero. 
In a probabilistic description this can be modeled by a stochastic force with zero mean. As the 
particle moves through the fluid, a hydrodynamic drag force sets in as a consequence of the random 
collisions with the fluid molecules that impedes the motion of the particle. Taking these two effects 
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into account, the equation of motion governing the particle is given by the Langevin equation 

=-??v(t)+f(t) . ( 1 ) 

Here v(t) and m are the instantaneous velocity and the mass of the particle, respectively, r/ is the 
friction coefficient, and f(t) represents the stochastic force. 

Before discussing the noise term further, let us stress that a stochastic description like equation 
( 1 ) assumes a separation of time scales between slow degrees of freedom (in this case the particle) 
and fast degrees of freedom (the fluid molecules). Usually the fast degrees of freedom represent the 
environment in which the slow degrees evolve. 

For the stochastic noise the following two assumptions are usually made: (1) the collisions of the 
fluid molecules with the particle are uncorrelated in time, and ( 2 ) the random forces are drawn 
from a Gaussian probability distribution with zero mean: 

(f(b) = 0, (2) 

-1') , (3) 

where v and ^ denote Cartesian coordinate indices and 7 is the strength of the noise. While these 
assumptions are perfectly sound for a Brownian particle, one should make sure for each situation 
whether they are indeed meaningful for the problem at hand. 

The noise term 7 , which measures the strength of the fluctuations, and the friction coefficient 7 , 
which describes dissipation, are related through the fluctuation-dissipation theorem (called Einstein 
relation in the present context): 


7 = 2/cbT?? , (4) 

where T is the temperature of the environment and /cs is Boltzmann’s constant. 

If the inertia time is larger than m/rj, inertial effects can be disregarded, and the particle’s 
motion is then governed by the overdamped Langevin equation 

= F(r(f))+f(t) , (5) 

where we have added a conservative force F(r) = — VU(r) originating from a potential U(r). 
Equation (5) is in fact the starting point for Brownian dynamics as well as for Langevin molecular 
dynamics used in the following. Note that imposing the fluctuation-dissipation theorem ensures that 
this stochastic differential equation correctly describes the equilibrium properties of our particle. 


2.2. Classical systems governed by a Hamiltonian 

The overdamped Langevin equation (5) is readily generalized to situations with N degrees of 
freedom (spatial coordinates of the particles) governed by a Hamiltonian iL[{r}] that is a functional 
of the position vectors r of the different slow degrees of freedom. In that case one has the following 
set of coupled stochastic differential equations 


dri{t) 


5H[ri{t)] 

5ri{t) 


+ fi(b ) 


( 6 ) 


where i = 1,... ,N labels the different particles. As before, the effects of the fast degrees of freedom 
coming from the environment are captured by the stochastic force fj(f) acting on particle i. These 
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stochastic forces have again zero mean, (fi(t)) = 0, and are taken from a Gaussian distribution, 
thereby fulfilling the Einstein relation (in d spatial dimensions) 

{^i{t) ■ fj(t')) = 2drjkBT5ijS{t - t') . (7) 

This set of equations guarantees that in equilibrium the system samples configurations with the cor¬ 
rect canonical equilibrium probability distribution, ~ ^-H/ksT^ However, out of equilibrium, which 
is our primary interest in the following, it is important to note that different microscopic dynamics 
might yield different results. This is even true in cases where the same probability distribution is 
ultimately obtained in equilibrium. A detailed comparison between different implementations of 
the system’s dynamics might therefore be in order if one wishes to ascertain that out-of-equilibrium 
results are not merely artefacts of the dynamics (see, e.g.. Ref. [6]). 


3. Langevin molecular dynamics for interacting vortex matter 


3.1. Elastic line model 

We consider in the following N interacting elastic lines in a disordered environment. Having in 
mind magnetic flux lines in high-Tc superconductors, we remark that this situation corresponds 
to the extreme London limit with the superconducting coherence length much smaller than the 
London penetration depth [7]. Vortex lines are then described by their trajectories rj(z), where 
z denotes the direction of the applied external magnetic field, and the two-dimensional vector rj 
indicates the xy position of line i at height z. As an immediate consequence of this description, 
magnetic flux lines in this model cannot form loops or overhangs, since rj(z) has to be unique. This 
is a reasonable assumption, as long as the vortex energy per length (i.e,. the elastic line tension) 
remains large compared to the energy scale of thermal fluctuations. 

The effective Hamiltonian of this system is written as a functional of the vortex line trajectories 
with an extent of L in the z direction and consists of three competing terms: the elastic line energy, 
the external potential due to disordered pinning sites, and the repulsive vortex-vortex interaction: 


HN[ri] 


_ N 




dri{z) 

dz 


2 N L 

dz + '^ VD{ri{z)) dz 
i=i do 



( 8 ) 


Consistent with the extreme London limit, the repulsive vortex-vortex pair interaction is set purely 
in-plane between different flux line elements. Below we will provide expressions for the attractive 
pinning potential and repulsive interaction appropriate for high-Tc superconductors. 

In order to employ a Langevin molecular dynamics algorithm to simulate the vortex line dy¬ 
namics, we discretize the system into layers along the z axis. Forces acting on the flux line vertices 
can then be derived from the properly discretized version of the Hamiltonian (8). We subsequently 
proceed to numerically solve the set of coupled overdamped Langevin equations 

ari(f,z) _ 5HN[ri{t,z)] , ^ 

- Sc.Xz) + (S) 

As discussed in the previous section, the fast, microscopic degrees of freedom of the surrounding 
medium are captured by thermal stochastic forces fj(t, z), modeled as uncorrelated Gaussian white 
noise with vanishing mean and fulfilling the Einstein relations (7) in d = 2 dimensions. This 
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guarantees that the system relaxes to thermal equilibrium with a canonical probability distribution 
in absence of an external driving force F^. In type-II superconductors such an external driving force 
stems from applied external currents via the Lorentz force. 

The time integration is performed via simple discretization of equation (9) [8]. In Ref. [6] we also 
considered the situation of small inertia times. The Langevin equation then acquires an additional 
inertial term and the time stepping algorithm becomes slightly more complicated [9]. 


3.2. Implementation for disordered type-II superconductors 

High-Tc superconducting materials are layered compounds and highly anisotropic: the lattice con¬ 
stant in the crystallographic c direction is much larger than the in-plane ones along the a and b 
directions; correspondingly, the effective charge carrier masses Mab and Me differ as well for the 
different directions. Henceforth we always assume that the magnetic field is aligned with the ma¬ 
terial’s crystallographic c direction, and the material properties discussed below are given for this 
configuration and assigned the in-plane index ab. When discretizing the system into layers along 
the z axis we choose the crystallographic c axis unit cell size cq as layer spacing [10, 11]. 

This anisotropy also shows up in the expressions for the model parameters when using the 
effective elastic line model (8). The elastic line stiffness or local tilt modulus is given by ei ~ 
r“^eo In {Xah/^ab)- Here Xab is the London penetration depth and ^ab is the coherence length, both 
in the ab plane, whereas T”^ = Mab/Me denotes the effective mass ratio. The energy per length eo 
is given by cq = {(/o/ATiXab)"^, with the magnetic ffux quantum (/q = /ic/2e. Finally the Bardeen- 
Stephen viscous drag coefficient y in the Langevin equation (9) \s y = c/q/27 :where pn 
represents the normal-state resistivity [12]. 

For our implementation of the Langevin molecular dynamics method for disordered type-II 
superconductors [6, 13-15] we also require expressions for the vortex-vortex interaction as well as 
for the attractive pinning potential. The repulsive in-plane vortex-vortex interaction is given by 

V{r) = 2eoKo{r/Xab) , (10) 


where Kq denotes the zeroth-order modified Bessel function. This is essentially a logarithmic re¬ 
pulsion that is exponentially screened at the scale Xab- In order to avoid artefacts due to periodic 
boundary conditions in directions perpendicular to z, we cut off the interaction at 5Aa;,. Pinning 
sites are modeled by randomly distributed smooth potential wells of the form 


Va{r,t) 



^a) 


1 — tanh 




( 11 ) 


where p > 0 is the pinning potential strength, and z^ and denote the z position and in-plane 
location of pinning site a, respectively. The pinning potential width bo will be our unit of length, 
whereas energies will be measured in units of eo^o- Th® full pinning potential VD{r,t) is then 
obtained by summing over all localized pinning sites in the system: 


No 

yD{r,t) = 14 (r,t) . 

a=l 


( 12 ) 


For the case of random point defects, which is the only case considered in this paper, these pinning 
centers are randomly distributed and chosen independently for each layer. If one were to instead 
consider columnar defects aligned parallel to the magnetic field along the z direction, then one 
would repeat the same spatial distribution pattern for each layer [6, 14, 16]. 
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The numbers we use for our model parameters are consistent with the actual values for YBCO 
[3]. We set the pinning potential width bo = 35A and also choose cq = 6o- The superconducting 
coherence length is then given by ^ab = 0.3 6o) whereas the London penetration depth is Xab = 346o- 
The effective mass anisotropy ratio is T”^ = 0.2 and the normal-state resistivity near Tc becomes 
Pn ~ 500 pQ cm. With this, the numerical values of all model parameters can be calculated. It 
then also follows that our fundamental temporal unit is to = rjbo/eo ~ 18 ps. In the following we 
shall measure all time intervals in units of to- 


3.3. Quantities of interest 

A large range of quantities of interest are at one’s disposal when studying the out-of-equilibrium 
properties of vortex matter. Some of these quantities are especially useful to characterize the steady- 
state properties in the presence of an external driving force, others allow us to capture the complex 
relaxation processes in situations where the system is initialized far from a steady state. 

A typical experimentally accessible quantity is the current-voltage (TV) characteristics. From 
a Langevin molecular dynamics study an equivalent description is obtained from a graph relating 
the driving force to the mean vortex velocity v. Indeed, from Faraday’s law the mean vortex 
velocity, which results from the velocities of each line element in direction of the driving force, is 
related to an induced electric field E = B x v/c and therefore to a voltage drop across the sample. 
A relation between the driving force F^ and an applied external current density j follows from the 
Lorentz force: F^ = j x 

The thermal spatial fluctuations along the vortex lines can be captured by the radius of gyration. 
Defined as Vg = \J{(vAz) — (rj))^), this quantity is in fact nothing else than the average root mean- 
square displacement from the mean lateral positions of the lines. Here, (• • •) indicates an average 
over all line elements of line i as well as an average over all N lines and different realizations of the 
disorder and the noise. 

A third way to explore steady-state properties of driven vortex matter is to study the pinning 
time statistics. This is done by calculating the distribution of flux line element dwelling times at 
defect sites. The dwelling time at pins is obtained by monitoring the positions of flux line elements 
and recording the instances where they enter and leave the attractive pinning centers. 

Non-equilibrium relaxation processes are best studied by preparing a system in a specific state 
before suddenly changing a system parameter [5]. This sudden change, often referred to as a quench, 
can be realized in many ways, depending on the physical system under investigation. For disordered 
vortex matter the following three types of quenches can be realized experimentally: temperature 
quenches, magnetic field (or vortex density) quenches, and external current (i.e., driving force) 
quenches. After such a quench the system ends up in a state far from stationarity. The resulting 
relaxation processes can then be monitored through various quantities. 

Many previous studies of glasses (structural glasses or spin glasses) as well as of systems with 
glassy-like dynamics have revealed that relaxation and related aging processes are best investigated 
through two-time quantities [5]. Indeed, as a system after the quench is far from stationarity, time 
translation invariance is broken and the properties of the system change with the time elapsed 
since the quench. As a result an observable that depends on two times s and t > s is not only a 
function of the time difference t — s but depends in more complicated ways on these two times. 

In the context of interacting vortices in a disordered environment an impressive list of two- 
time quantities has been discussed [6, 14, 15, 17-21]. Some of these quantities contain information 
about local thermal fluctuations, whereas others provide insight into the time evolution of the 
global structure of the flux line configuration. We refrain from discussing all these quantities, but 
instead only introduce those used and displayed in the following concise overview: The height- 
height autocorrelation function (sometimes simply called the roughness) probes the local thermal 


6 



November 5, 2015 Molecular Simulation langevin 


fluctuations of vortices about their mean lateral position. It is defined as 

C'(t,s) = ((ri(t, 2 ;) - (rj(t)))(rj(s, 2 ;) - (rj(s)))) , (13) 

where the averages are again taken over the flux line elements of all lines as well as noise and 
disorder realizations. The average square distance between a vortex line element’s position at the 
two times s and t is measured by the two-time mean square displacement 

B{t,s) = {{ri{t,z)-ri{s,z))‘^) . (14) 

This quantity contains information on the decay and formation of global structures. 


4. Disordered vortex matter out of equilibrium 

We illustrate the rich physics of interacting magnetic vortices in disordered type-II superconductors 
that is governed by competing energy and hence length and time scales, through various character¬ 
istic examples that address both fluctuations in non-equilibrium steady states [6, 13, 15] as well as 
complex out-of-equilibrium relaxation scenarios initiated from very different starting configurations 
[6, 14, 15]. In this brief overview, we focus on systems with comparatively weak magnetic fields 
and hence low vortex density, and at low temperatures. Consequently the equilibrium stationary 
configurations in the presence of attractive point-like pinning centers likely reside in the disorder- 
dominated Bragg glass thermodynamic phase, wherein spatial positional order is disrupted at long 
length scales. At the cost of elastic energy, the flux lines attempt to accommodate as many localized 
defect sites as possible, which renders their trajectories through the sample quite rough [22]. 


4.1. Driven vortex lines in disordered environments 

When sufficiently large external currents are applied, the vortices detach from the pinning centers 
and freely flow through the sample. They then tend to reorder into regular triangular moving 
lattices. At zero temperature, the transition from the pinned, glassy phase with drift velocity 
u = 0 to the moving state with u > 0 is sharp, and represents a continuous non-equilibrium phase 
transition [23]. At finite temperatures, this transition becomes smoothened out as the external 
driving force is increased. This is shown in Fig. 1 (top panel) for a Langevin molecular dynamics 
simulations of a system of V = 16 vortex lines (discretized into L = 100 6o layers transverse to the 
magnetic field) driven in a direction perpendicular to the external magnetic field, at temperature 
T = 0.002 eofeo/^B Ri 10 K and with pinning strength p = 0.05eo- Below the depinning threshold 
Fc ~ 0.002 eo/6o) the flux lines hardly move, but display thermal fluctuations; the system remains 
superconducting. For larger drives, the mean vortex speed scales linearly with the driving force. 
This indicates Ohmic dissipation, since F^ oc j and v oc E, the induced electric field. Thus the 
freely flowing vortex system represents a normal-conducting phase. 

The Langevin molecular dynamics simulations allow tracking of each individual vortex line, and 
hence a detailed analysis of the statistics of the transverse line fluctuations that are induced both 
thermally and through interactions with localized pinning sites. At large driving forces F^ S> Fc, 
the attractions exerted by point defects become irrelevant, and thermal energies are also minute 
compared to the work associated with the drive. Consequently, the flux lines become straightened 
out, as is clearly seen in Fig. 1 (bottom panel), which displays the mean vortex radius of gyration 
as function of F^- In the disorder-dominated Bragg glass phase, the gyration radius is enhanced by 
about a factor of three, indicating the roughness of the vortex trajectories. As Fa approaches Fc, 
the drive pulls more line segments away from the attractive defects, which increases the radius of 
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Figure 1. Mean vortex velocity v (top), in units bo/to, and flux line gyration radius Vg (bottom), in units of bo, in the non¬ 
equilibrium steady state vs. driving force in units of eo/6o, for a system of = 16 interacting vortices of length L = 100 6o 
at temperature T = 0.002 epbo/fes and for point pinning center strength p = 0.05 eq. The top graph directly maps onto the 
current-voltage characteristics for superconductors; the location of the gyration radius maximum in the bottom graph indicates 
the depinning threshold [15]. 


gyration. A distinct maximum of Vg is observed right at the depinning threshold, a clear signature of 
the aforementioned zero-temperature non-equilibrium phase transition, for which rg{Ffi —>■ Fc) 
oo would display a critical divergence [ 6 , 15]. 

Even more detailed information on the pinning features of flux lines to variable point disorder is 
offered by careful analysis of the dwelling time statistics. Over many simulation runs, we recorded 
the durations r that each line element of a single driven vortex (with L/ 6 o = 200 segments) spent 
within a lateral distance bo of its previous position, and gathered the resulting histogram P{t) 
[13]. We kept the driving force Fa = 0.002 eo/ 6 o fixed, but varied the temperature T and the 
width (standard deviation) w of an assumed Gaussian distribution for the pinning strengths p 
about the mean p. Since the point defects are correlated only over a distance bo in our system, 
the collective Larkin-Ovchinnikov pinning length scale is ~ boiei/w)"^^^. The associated tpyical 
energy barrier then becomes ~ eib^/L^, from which one estimates the critical depinning force 
as Fc ~ El/{ bo Ll) k, rc‘^/^/(6o G^^) [24]. Employing extreme-event statistics arguments, Vinokur, 
Marchetti, and Chen argued that in the pinned phase, the pinning time distribution should obey 
a power law for large dwelling times r, P{t) ~ with an effective temperature-dependent 

scaling exponent a{T) ~ ksT/Ei [24]. 

This power-law behavior is indeed borne out by the Langevin molecular dynamics simulation 
data depicted in Fig. 2 (left panel, a) at large disorder variance w = 1.0 cq. In contrast, in the weak 
pinning regime {w = 0.001 cq), the vortices remain attached to the defects only for short times 
r < 1000fo) before they detach and freely flow through the system. The marked maximum of P{t) 
observed in the corresponding graph may serve to define a typical pinning time r ~ 700to- The 
dynamical phase diagram associated with the transition between the freely flowing and the pinned 
regime has been studied by Krauth et al by numerically extracting the structure factor exponent 
of an elastic string in the T —>• 0 limit [25] and in the finite temperature case [26]. The right panel 
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Figure 2. Left panel (a): Distribution of pinning times P{t) for single vortices with L/6o = 200 line elements moving through a 
system with N^/L = 100 randomly placed isolated point defects, whose pinning strengths are Gaussian distributed with mean 
/i = 0 and variance driving force = 0.002 temperature T = 0.002 cq In th® flowing phase (w = 0.001 cq), 

a characteristic dwell time can be extracted from the observed maximum. In the pinned phase {w = 1.0 eo), F(t) decays 
algebraically with r. Right panel: (b) Temperature dependence of the power-law scaling and effective decay exponent o:(T) for 
the pinning time distribution P{r) (with T measured in units of eo bQ/ks)- (c) Similar results, but obtained for a continuous 
Gaussian disorder landscape correlated over distance ^ = 6o with otherwise identical parameter values [13]. 


(b) of the figure explores the temperature dependence of the effective decay exponent a{T)] as the 
inset for the low-temperature regime T < 10“^ f-obo/ks confirms, one indeed finds a roughly linear 
dependence a(T) ~ 0.0417 -I-3340T, albeit with an apparent non-zero intercept as T —>• 0 [13] that 
is not captured by the scaling theory. We also investigated the dwelling time distribution where 
the discrete point defects are replaced by a smooth, continuous disorder landscape, in order to 
make more direct contact with the theoretical setting in Ref. [24]. To this end, we drew the pinning 
potential strength p at each node of a square lattice with spacing bo from a Gaussian distribution 
with vanishing mean at variance w'^, and constructed the continuous landscape from the resulting 
grid points through straightforward bi-linear interpolation. In the pinned state, one again observes 
algebraic decay of P{t), but actually less cleanly, as shown in the the right panel (c) of Fig. 2; 
also, the typical dwelling time in this situation turns out to be r ~ 100 to > indicating a different 
overall energy scale as compared with the discrete pin case. At low temperatures, we found the 
linear relationship a{T) « 0.2542 -|- 2019 T (inset) [13]. 


4.2. Relaxation dynamics of vortex lines in disordered type-II superconductors 

Next we proceed to investigate the non-equilibrium relaxation kinetics of vortex matter in disor¬ 
dered type-II superconductors starting from various initial configurations that are quite distinct 
from the stationary states that may be reached after long time durations. In this brief overview, 
we report Langevin molecular dynamics simulation data for mutually interacting flux lines sub¬ 
ject to randomly distributed attractive point pinning centers, and probe the relaxation kinetics 
through the various two-time quantities introduced in section 3.3. We first study simulations for 
systems where initially perfectly straight vortex lines were placed at random positions in the 
sample, and then allowed to relax towards Bragg glass equilibrium states at low temperatures 
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t/s 

Figure 3. Relaxation of the two-time mean-square displacement s) for flux lines of length L = 640 bo subject to attractive 
point defects of strength p = 0.05 eo? at temperature T = 0.002 eo5o/^s? and for various waiting times s as indicated, starting 
from initially straight vortices placed at random locations: (a) in the absence of mutual interactions; (b) in the presence of 
repulsive forces. Simple aging scaling collapse of the data can only be achieved for small waiting times s, with the aging 
exponents b 0.725 and b ^ 1.54. [6]. 


T = 0.002 eobo/ks ~ lOK at vanishing driving force = 0. We shall attempt to fit the data to the 
simple aging scaling form A{t, s) ~ FA{t/s) for two-time observable A with the aging exponent b. 
Indeed, in the absence of disorder and mutual vortex interactions, this setup leads to simple aging 
scaling for both the two-time mean-square displacement B(t,s) and the two-time height-height 
autocorrelation function C{t,s) with the Edwards-Wilkinson scaling exponent b = 0.5 [6, 19, 27]. 

In the presence of attractive point-like pinning centers, however, such simple aging scaling behav¬ 
ior is not observed, although time translation invariance is manifestly broken. Figure 3 displays data 
for the two-time mean-square displacement for various waiting times s, both for non-interacting 
flux lines (a) and mutually repelling vortices (b) with L/bo = 640 segments. Simple aging data 
collapse can only be obtained at short times s, with the scaling exponents b ~ 0.725 and b ~ 1.54, 
respectively. Non-interacting vortices strive to optimize the balance between pinning energy gains 
and elastic stretching energy losses. Interacting flux lines in addition become ‘caged’ through the 
repulsive interactions with their neighbors, which drastically reduces lateral line fluctuations and 
displacements. These same competitive effects are visible also in the simulation data for the height- 
height autocorrelation function C{t, s). Simple aging scaling again ensues only for very short wait¬ 
ing times s, with, e.g., a scaling exponent b ~ 0.716 for non-interacting vortices as indicated in 
Fig. 4(a). The competition between elastic forces and pinning induce remarkably complex and 
even non-monotonic relaxation features. When the repulsive in-plane pair interactions between 
flux line elements are taken into account, the two-time height-height autocorrelation data display 
characteristic two-step relaxation behavior, see Fig. 4(b), as is also frequently seen in structural 
and spin glasses [28]. In the ‘/3’ relaxation time window, C{t, s) changes hardly at all, and the data 
appear to satisfy time translation invariance; in the ultimate ‘a’ relaxation regime, the normalized 
height-height autocorrelation function finally decays, albeit very slowly [6]. 
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Figure 4. Relaxation of the two-time height-height autocorrelation function (^(t, s) for flux lines of length L = 640 6o subject 
to attractive point defects of strength p = 0.05eo, at temperature T = 0.002 eo^o/^s, and for various waiting times s as 
indicated, starting from initially straight vortices placed at random locations: (a) in the absence of mutual interactions; (b) in 
the presence of repulsive forces. [6]. 


We pause for a moment to mention the work of Bustingorry et al [17, 18] in which relaxation 
processes in a related system were studied using Langevin molecular dynamics. In that work both 
attractive and repulsive pinning sites were present in the system. This modification yields important 
changes in the relaxation processes. For instance, instead of observing the two-step glassy-like 
relaxation shown in Fig. 4, a simple aging scaling of the autocorrelation prevails. 

Experimentally, the initial configuration of randomly placed straight vortices studied above 
cannot be realized. In order to address more realistic initial conditions, we therefore ran a series 
of numerical experiments wherein the system is allowed to relax for 10^ time steps, until at t = 
100,000 to a sudden change of an external experimental control parameter is implemented. We 
then monitored several observables to characterize the system’s non-equilibrium relaxation features 
following such quenches. A first example is depicted in Fig. 5: At the quench point, the system’s 
temperature was instantaneously raised from T = 0.002 eobo/ks to T = 0.008 eo^o/^B- Both the 
radius of gyration rg{t) and the fraction of pinned flux line elements ^{t), defined as the fraction 
of segments that reside within a distance 6o of the attractive point pins at time t, are seen to relax 
essentially exponentially towards the system’s new steady state at elevated temperature. For the 
gyration radius, the relaxation time in Fig. 5(a) is measured to be r ~ 34600 to- Asymptotically, the 
temporal evolution approaches that of an unperturbed vortex sample at the higher temperature. 
Careful inspection of Fig. 5(a) reveals a tiny dip immediately following the quench before the 
expected larger gyration radius value at elevated temperature is reached. Indeed, as confirmed 
in panel (b), at first a few line elements become thermally depinned by the sudden temperature 
increase, which initially slightly relaxes the disorder-induced vortex line roughening in the Bragg 
glass phase. At any rate, the exponentially fast relaxation prevents the emergence of a sizeable 
aging time window [14]. 

Considerably richer and more complex relaxation features are observed following a magnetic 
field quench at time r = 100,000 to j when the number of vortices in the system is either suddenly 
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Figure 5. Time evolution of (a) the radius of gyration rg{t), and (b) the fraction of pinned vortex line elements following a 
sudden temperature quench from T = 0.002 eo^o/Zc^ to T = 0.008 eo^o/Zcg at t = 100, 000 to for a system of = 16 interacting 
flux lines of length L = 640 5o subject to attractive point defects of strength p = 0.05 cq. For comparison, the time tracks at 
fixed initial and flnal temperatures are shown as well. The system relaxes exponentially fast to its new steady state [14]. 




Figure 6. Relaxation of the normalized two-time height-height autocorrelation function for various waiting times s (as indi¬ 
cated) at temperature T = 0.002eo5o/^s for flux lines of length L = 6406o subject to attractive point defects of strength 
p = 0.05 eo (a) at fixed magnetic field or flux density, with N = IQ flux lines in our simulation sample; (b) following a field 
down quench, suddenly reducing the vortex number from = 21 to 16; (c) following a field up quench, where 5 additional flux 
lines are added to the original N = IQ vortices; (d) as in (c), but displaying the data only for the 5 newly added flux lines [14]. 


reduced (down quench) or increased (up quench), in both cases maintaining a fixed temperature 
T = 0.002 eobo/ks (and a vanishing driving force = 0). We introduce the total elapsed simulation 
times r = r +1 and cr = r + s; since r will be much larger than the times s and t, both measured 
after the quench at r, the shift by the duration r needs to be taken into account in attempts to 
achieve simple aging scaling. Figure 6(a) displays the relaxation curves for the normalized two-time 
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Figure 7. Time evolution of the two-time height-height autocorrelation function for various waiting times s (as indicated) 
at temperature T = 0.002 eo^o/^s for driven flux lines of length L = 640 6o subject to attractive point defects of strength 
p = 0.05 eo, as the driving force is instantaneously reduced from = 0.015 eo/5o (i.e., from a moving non-equilibrium steady 
state) into the pinned Bragg glass phase at Fd = 0.001 eo/^0- The inset shows the data plotted vs. t — s, demonstrating the 
breaking of time translation invariance; the main panel indicates the best attempt to fit the data to a simple aging scaling form 
with b = 0.392, which however only partially works in a very limited time window [15]. 


height-height autocorrelation function with the magnetic field or vortex density held constant for 
various waiting times s = cr — r; the data essentially correspond to the longest waiting time curve 
in Fig. 4(b) and demonstrate that for the large duration r = 10® to considered here, the vortex 
system has basically reached the stationary regime for which time translation invariance holds: the 
data collapse when plotted against the time difference F — cj = t — s. In contrast, clear evidence 
of broken time translation invariance is visible in Fig. 6(b), where at r = 100, 000 to five randomly 
selected vortices were instantaneously removed from the system. One still observes the a-/3 two-step 
relaxation scenario, yet the system has not reached stationarity yet, akin to Fig. 4(b) [14]. 

Qualitatively very similar relaxation features are observed in the system of originally N = 16 
vortices, when at r = 100,000 to five new and initially straight flux lines are introduced at random 
positions in the sample, as seen in Fig. 6(c). Time translation invariance is again manifestly broken 
in the late a relaxation regime. By separately assessing the relaxation data for the two-time height- 
height autocorrelation function of just the newly added flux lines, depicted in Fig. 6(d), it becomes 
apparent that the distinction between Figs. 6(b) and (c) are caused by the marked non-monotonic 
behavior of the lateral line fluctuation relaxations of the added line population, which typically 
first show a strong increase followed by a much slower final decay. Only at very long waiting 
times s > 2^^ to do we observe the monotonic a-/3 relaxation features. When the new vortices are 
inserted into the sample, they experience strong repulsive forces from the originally present flux 
lines, which increases their range and thus facilitates their pinning to point defects. This in turn 
enhances transverse line fluctuations, until at ln(t — s) ~ 7 optimal pinning configurations have 
been reached, whence the final slow decay of the autocorrelation function ensues [14]. 

As a final set of numerical experiments, we initialize the vortex system with = 16 lines at 
temperature T = 0.002 eobo/ks and at sufficiently large external driving force = 0.015 eo/bo 
to ensure it resides well in the moving non-equilibrium steady state, see Fig. 1. After relaxing 
the system for a duration of 200,000 toj the drive is suddenly lowered to = 0.001 eo/bo, which 
forces the flux lines into the pinned, disorder-dominated Bragg glass phase. Subsequently, we again 
monitor various observables; as characteristic example, the relaxation data for the two-time height- 
height autocorrelation function C(t, s) are shown in Fig. 7. As in the previous magnetic field quench 
scenarios and for randomized initial conditions, simple aging scaling does not ensue for this quantity, 
but holds only approximately in a very limited time regime (with scaling exponent b ~ 0.392). It 
is worth mentioning that upon starting in the pinned glassy state, and suddenly increasing the 
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driving force into the moving state, the relaxation processes occur exponentially fast and quickly 
reach stationarity, with relaxation times of the order of 1200 fo [15]. 


5. Conclusion 

Langevin molecular dynamics is a powerful numerical method for the investigation of interacting 
systems with many degrees of freedom and clear separation of time scales. In this paper we have 
illustrated this through the study of interacting vortex matter in different far-from-equilibrium 
situations. We focused on values for the system parameters that are realistic for high-Tc super¬ 
conductors. When driving the vortices through the sample by means of an external current, non¬ 
equilibrium steady states and a (zero temperature) non-equilibrium phase transition emerge. As 
the Langevin molecular dynamics simulations allow to track each individual vortex line, a de¬ 
tailed investigation of the steady-state properties is possible through the measurement of a variety 
of quantities. A second experimentally relevant situation arises when a system is brought out of 
equilibrium through a rapid change (quench) of some of the external system control parameters. 
Vortex matter in a disordered environment provides an extremely rich and technologically impor¬ 
tant model system, with a variety of such perturbations that are experimentally feasible, namely 
temperature, magnetic field, as well as current quenches. Detailed investigations of these different 
situations allow us to gain a rather complete understanding of the relaxation processes that take 
place when interacting vortex matter is forced out of thermal equilibrium. 
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